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Abstract 

We give an algorithm, with R code, to minimize the multidimensional scaling stress 
loss function under the condition that some or all of the fitted distances are smaller 
than given upper bounds. 
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1 Introduction 

In this paper we study the minimization of stress, defined as 

inn 

stress(X) := - - (UfiX)) 2 

4 i=1 j =l 

over n x p configuration matrices X. Here W = {wy,} and A = {hp} are given symmetric, 
hollow, and non-negative matricesof, respectively, weights and dissimilarities. The d l3 {X) are 
Euclidean distances between the points with coordinates in the rows of X. Thus 

dij( X) := (xi - Xj)'(xi - xj) = (e* - efi'XX '(e* - efi = tr X'AijX. 

Here e* and e 3 are unit vectors, with all elements equal to zero except for the single element i 
or j, which is equal to one. Also A VJ := (e; — efiici — efi)'. 

The problem of minimizing stress is well-studied (see, for example, Borg and Groenen (2005)) 
. In this paper we analyze the problem of minimizing stress over all configurations X for 
which dij(X) < aij for some or all pairs i,j , where the are positive bounds. As a special 
case we have MDS from below , where we require difiX) < 5ij for all i,j. 

Note that constraints of this type are not exactly independent. If we require d l3 (X) < a t] and 
dik{X) < an- then, by the triangle inequality, we automatically require d 3 k{X) < a t] + cq fc . 
On the other hand the constraints are always consistent, because they will be satisfied if we 
take X small enough. 

The constraints difiX) < « y -, or equivalently rf|-(X) < a 2 -, define a convex set in configuration 
space, but stress is not convex. We can, however, use basic smacof theory (De Leeuw (1977)) to 
majorize stress by a convex quadratic. This can be used to define a simple iterative algorithm. 
In each iteration we majorize stress, and then minimize the convex quadratic majorizer over 
configurations satisfying the convex quadratic constraints dfj(X) < afj. Compute a new 
majorization at the minimizer of the majorized problem, and so on. These are the two steps 
in the majorization or MM approach to optimization (De Leeuw (1994), Lange (2016)). To 
minimize a convex quadratic under convex quadratic constraints we use the dual Newton 
method from De Leeuw (2017). 

1.1 Simplifying stress 

A simplified expression for stress, first used in De Leeuw (1993), is 

1 K r - 

stress(x) = - J2 w k( 4 - sJx'Akx) 2 . 

k=l 

We have put the elements below the diagonal of W and A in a vector of length \n{n — 1). 
Also x := vec(A) and the Ak are the direct sums of p copies of the corresponding A tJ . Now 
expand the square, and assume without loss of generality that \ J2k=i w kd 2 = 1. 
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Then 


stress(x) 


with V : = Ef=i w k A k . 


1 - 


J2w k 8 k y/x'A k x + 

k =1 



Suppose V = KA 2 K’ is a complete eigen-decomposition of V. Some care is needed to handle 
the fact that V is singular, but under the assumption of irreducibility (De Leeuw (1977)) 
this is easily taken care of. Change variables with z = AzK'x. Then 


stress^) 



1 

2 


/ 

z z , 


with A k = A 2 K'A k KA a, so that Y^ k =i w kA k = /. 


1.2 Smacof Notation and Theory 

We have see in the previous section that we can write 

stress(x) = 1 — p(x) + -x'x, 


where 


I\ 

p(x) ■■= J2 w ^kd k (x), 

k =1 


with distances d k (x) := y/x'A k x, and the A k positive senri-definite matrices that add up to 
the identity. By Cauchy-Schwarz, 


p(x) = x'B(x)x > x'B(y)y, 


where 

B(x) := £ l Wk m At - 

If we define the Guttman transform if x as T(x) := B(x)x, then for all x 

(. stress)(x ) = 1 — x'T(x) + ^ x'x = (1 — ^T(x) / r(x)) + ^(x — T(a;)) / (a; — T(x)), 
and for all x, y 

(stress) < 1 - x'Y(y) + ^ x'x = (1 - ^T(j/)T(y)) + ^(x - Y(y)'(x - T(y)), 

In this notation a smacof iteration is x^ k+1) = T^*-^). 
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1.3 Necessary Conditions 

We look at the necessary conditions for an optimum of the bounded MDS problem. The 
Lagrangian is 

1 K 

C(x, A) := stress(x) + - Y A k (x'A k x - S 2 k ) 

^ k =1 

If x is feasible and optimal for the problem of minimizing stress from below then there is a 
A > 0, not all zero, such that 

V x L{x, A) = | f 1 + A k A^j x - T(x) = 0 

and the multipliers A and constraint functions are complementary, i.e. A k {x'A k x — 8 k ) — 0 
for all k. Compare this with the necessary condition for a local minimum in regular smacof, 
which is simply x = 

Gamma{x). 

It should be noted that x = T(x) implies p(x) = x'x, and consequently at a stationary 
point cr(x) = 1 — \x'x , which implies x'x < 2. Thus Y k =i w kd\{x) < 2, which implies 

dk(x) < yj2/wk■ Even in regular smacof, without bound constraints, the distances at a 
stationary point are bounded, and imposing the constraints x'A k x <2 /w k still allows the 
unconstrained smacof solution. If the minimum stress is non-zero, then none of the constraints 
will be active, and all the distances will be strictly smaller than the a priori bounds. 


2 MDS from Below 


Classical MDS (Torgerson (1958)) can be interpreted as minimizing the loss function strain , 
defined as 

strain(X) := ^tr J n ( A 2 - D 2 (X))J n ( A 2 - D 2 (X)), 

over all n x p centered configuration matrices X , where J n is the centering operator J n = 
I n — -ee', A 2 is a matrix with squared dissimilarities, and D 2 (X) is a matrix of squared 
Euclidean distances. 


The loss function interpretation is based on the identity 

-ij n D 2 (X)J n = AX', 


which implies 
with C := — \J n X 2 J n . 


strain(A") = tr (C — XX') 2 , 


If C has signature (n + ,no,n_) then we can write C = C — C_ with both C and C_ positive 
semi-definite, of ranks n + and n_, and with tr CC_ = 0. C is the projection of C on the 
convex cone of positive semi-definite matrices, and C = —{C — C) is the negative of the 
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projection on its polar cone. In terms of the eigenvalues and eigenvectors of C we can write 
C = KAK + KAK . with overlining used for the positive eigenvalues and their eigenvectors, 
and underlining for the negative ones. Thus 


strain(X) = tr (C-C- XX') 2 = tr (C - XX') 2 + tr C 2 + 2tr X'CX, 


and 

min strain(X) > mintr (C — ZZ') 2 + tr C 2 + 2 rnintr Y'CY = rnintr (C — ZZ') 2 + tr C 2 . 

Thus if p > n + we choose Z so that ZZ' — C, if p < n + we use the p largest eigenvalues and 
corresponding eigenvectors of C to find ZZ'. 

Now consider the squared distances d 2 j(C) = Ca + Cjj — 2and the squared distances 
dpX r ) = V_, - Xjs) 2 - Then for all i,j 

4(AT) < d%( X 2 ) <•••< d%(X n+ ) = 4(C). 

Thus solutions are nested and we use quadratic approximation from below of the positive 
semi-definite part of C (De Leeuw and Meulman (1986)) 

The idea of approximation from below, which is natural in the case of classical MDS, can also 
be applied to the minimization of stress, using the coinstraints dij(X) < Sij for all % and j. 

This may seem somewhat perverse, because why would we ever impose such constraints ? In 
this context, however, it makes sense to interpret the S t] as bounds instead of dissimilarities. 
The problem is to locate objects in space in such a way that their distances do not exceed 
their upper bounds, and are as close to these bounds as possible. Thus the upper bound 
problem is perhaps more like a location problem than a multidimensional scaling problem. 


3 Examples 

3.1 Equidistances 

Our first example has n = 10 with all S t] equal. The optimal smacof solution in two 
dimensions needs 131 iterations to arrive at stress 0.1098799783. Approximation from below 
uses 30 iterations and finds stress 0.1520077612. The two configurations are in figure 1. 
Note they produce two different local minima, one with nine points equally spaced on the 
circle and a tenth point in the center, and the other with ten points equally spaced on the 
circle. At the solution of the bounded problem there are 10 active constraints, which have 
dij(X) = SijiX) = 1 . 


5 




CD 

d 



8 

1 



d 

7 

10 

5 




3 




CM 





C\J 




6 


o 

4 


2 

C\J 

d 





CM 

o 




E 



4 

9 


E 

o 




■a 

CM 





T3 


q 




d 

1 




2 



o 


6 


CD 


7 




0.4 

J_ 

8 


9 


d 

1 



10 5 



i 


1 





CD 

— d 

i 

CM 

d 

CM 

— d 

i 

1 

0.6 



1 1 

-0.4 

o 

d 

1 1 

0.2 0.4 


dim 1 dim 1 


Figure 1: Equidistance Data, Regular Left, From Below Right 




dissimilarities dissimilarities 

Figure 2: Equidistance Data, Regular Left, From Below Right 

To illustrate complementary slackness we give the values of the constraint functions d'j 3 (X)— <5 7, 
which are all non-positive, and the values of the Lagrange multipliers, which are all non¬ 
negative. We see that a constraint function that is equal to zero means a positive Lagrange 
multiplier, and a constraint function that is negative means a zero Lagrange multiplier. 

mprint (h2b$constraints, d = 6, f = "+") 


## [ 1 ] 
## [ 8 ] 
## [15] 
## [ 22 ] 
## [29] 
## [36] 


-0.172287 -0.309015 
-0.451058 +0.000000 
-0.093098 -0.344574 
-0.451053 -0.172291 
-0.093098 -0.344576 
+0.000000 -0.172292 


-0.172287 +0.000000 
+0.000000 -0.093093 
-0.344576 -0.451059 
-0.172292 -0.172283 
-0.309023 -0.309009 
-0.451053 -0.172292 


-0.309015 +0.000000 
-0.451058 -0.451058 
+ 0.000000 + 0.000000 
+0.000000 -0.451058 
+0.000000 -0.172288 
-0.172287 +0.000000 


-0.451059 

-0.172283 

-0.309023 

-0.344574 

-0.451054 

-0.451055 
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## [43] -0.344584 -0.093096 -0.093096 
mprint (h2b$multipliers, d = 6, f = "+") 


## [ 1 ] 
## [ 8 ] 
## [15] 
## [ 22 ] 
## [29] 
## [36] 
## [43] 


+ 0.000000 

+ 0.000000 

+ 0.000000 

+ 0.000000 

+ 0.000000 

+0.088508 

+ 0.000000 


+ 0.000000 

+1.220853 

+ 0.000000 

+ 0.000000 

+ 0.000000 

+ 0.000000 

+ 0.000000 


+ 0.000000 

+1.220844 

+ 0.000000 

+ 0.000000 

+ 0.000000 

+ 0.000000 

+ 0.000000 


+0.088502 

+ 0.000000 

+ 0.000000 

+ 0.000000 

+ 0.000000 

+ 0.000000 


+ 0.000000 

+ 0.000000 

+0.088508 

+1.220844 

+1.220840 

+ 0.000000 


+0.088503 

+ 0.000000 

+0.088477 

+ 0.000000 

+ 0.000000 

+1.220840 


+ 0.000000 

+ 0.000000 

+ 0.000000 

+ 0.000000 

+ 0.000000 

+ 0.000000 


3.2 Dutch Political Parties 1967 

As the next illustration we use data from De Gruijter (1967), with average dissimilarity 
judgments between Dutch political parties in 1967. The data are 

## KVP PvdA VVD ARP CHU CPN PSP BP 

## PvdA 5.63 

## VVD 5.27 6.72 

## ARP 4.60 5.64 5.46 

## CHU 4.80 6.22 4.97 3.20 

## CPN 7.54 5.12 8.13 7.84 7.80 

## PSP 6.73 4.59 7.55 6.73 7.08 4.08 

## BP 7.18 7.22 6.90 7.28 6.96 6.34 6.88 

## D66 6.17 5.47 4.67 6.13 6.04 7.42 6.36 7.36 

The optimal smacof solution in two dimensions needs 319 iterations to arrive at stress 
0.044603386. Approximation from below uses 32 iterations and finds stress 0.0752702106. 
At the solution there are 10 active constraints, but as the Shepard plots in figure 4 show, 
the active constraints now do not correspond with the largest dissimilarities, and the two 
configurations in figure 3 are quite different. Although that could, of course, also be because 
they correspond with two different local minima, which mostly differ in the location of the 
D66 party (a neoliberal party, full of civilized protest). 
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Figure 3: De Gruijter Data, Regular Left, From Below Right 




Figure 4: De Gruijter Data, Regular Left, From Below Right 

In the next analysis we require that all distances are less than or equal to 8.13, the maximum 
dissimilarity. In order words, the maximum distances must be less than or equal to the 
maximum dissimilarity. 

The optimal solution under these constraints in two dimensions needs 70 iterations to arrive 
at stress 0.0537645693. At the solution there are 6 active constraints, i.e. six distances equal 
to 8.13. 
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Figure 5: De Gruijter Data, Distances Bounded by Maximum Dissimilarity 

And finally an analysis in which we require that distances between ARP, CHU, KVP (the 
Christian Democrat parties) are less than 0.1, and the distances between PvdA, PSP, CPN 
(the leftist parties) are also less than 0.1. This will effectively collapse them in the configuration 
plot. 

The optimal solution under these constraints in two dimensions needs 183 iterations to arrive 
at stress 0.0494598864. At the solution all 1 constraints are active. 
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Figure 6: De Gruijter Data, Imposing Small Distances 


3.3 Ekman Color Data 

The final example are the color data from Ekman (1954). The optimal smacof solution in 
two dimensions needs 26 iterations to arrive at stress 0.0172132469. Approximation from 
below uses 16 iterations and finds stress 0.0274106473. At the solution there are 15 active 
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constraints. Again the bounded solution is a shrunken version of the regular smacof solution, 
and figure 8 shows that the active constraints correspond with the largest dissimilarities and 
distances. 
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Figure 7: Ekman Data, Regular Left, From Below Right 
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Figure 8: Ekman Data, Regular Left, From Below Right 


Figure 9 shows that the two configurations in figure 7 may look to be proportional, they 
really are not. There are many deviations from proportionality if we compare the fitted 
distances in the regular and bounded smacof solutions. 
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4 Discussion 

If p — n — 1 our multidimensional scaling problem is full-dimensional scaling or FDS (De 
Leeuw, Groenen, and Mair (2016)). In that case stress can be parametrized by using the 
scalar products C and we must minimize 

n n 

stress(C) := - \/ tr A ij C Y 

i= 1 j =1 

over C > 0 and tr AijC < ctij. This means that FDS with upper bounds means minimizing 
a convex function over a convex set, in fact solving a semidefinite programming problem. It 
also means that the concept of Gower rank starts playing a role (De Leeuw (2016b)). 

MDS from below can also be used with different loss functions. Minimizing 

n n 

\hj-dij( x )\ 

i =1 3=1 

over dij(X) < Sij is the same as maximizing the weighted sum of the distances, a convex 
function. Thus we can use tangential majorization (De Leeuw (2016a)), which means the 
majorization function is linear and has no quadratic term. 

Alternatively, the problem of minimizing 

n n 

£!><#S-4-P0l 

*=1 3=1 
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over dij(X ) < S t] means maximizing a convex quadratic form. Although the objective function 
is quadratic, we still need tangential majorization to get a convex objective to minimize. 

We do not address the problem in this paper of minimizing stress under constraints of the 
form (3ij < dij(X) < oiij. This is a fundamentally more complicated problem, because the 
lower bounds make the solution set non-convex. In fact, using lower bounds may make the 
solution set empty, because there may not be p-dimensional Euclidean distances satisfying 
the constraints. 

It will be of some interest to apply our algorithm to unfolding examples, where we could 
have upper bound constraints on the row distances, the column distances, or the row-column 
distances. Of course in the unfolding case lower bound constraints are at least as interesting, 
because they can in principle be used to avoid the trivial solutions that are so common in 
unfolding, but we have argued that computationally lower bound constraints are far more 
difficult to handle. 


5 Appendix: Code 

5.1 below.R 


suppressPackageStartupMessages 

suppressPackageStartupMessages 

suppressPackageStartupMessages 

suppressPackageStartupMessages 

source (" auxilary.R") 

source ("mdsUtils.R") 

source ( "nnnewton.R" ) 

source ( "qpqc.R" ) 

source ("smacof.R") 

source ( " smacofBelow.R" ) 


(library (mgcv, quietly = TRUE)) 
(library (MASS, quietly = TRUE)) 
(library (lsei, quietly = TRUE)) 
(library (numDeriv, quietly = TRUE)) 


5.2 auxilary.R 


mprint <- function (x, 

d = 6, 
w = 8, 
f = "") { 

print (noquote (formate ( 

x, 

di = d, 
wi = w, 
fo = "f", 
flag = f 
))) 
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> 


directSum <- function (x) { 
m <- length (x) 
nr <- sum (sapply (x, nrow)) 
nc <- sum (sapply (x, ncol)) 
z <- matrix (0, nr, nc) 
kr <- 0 
kc <- 0 

for (i in l:m) { 
ir <- nrow (x [ [i]]) 
ic <- ncol (x[ [i]] ) 

z [kr + (l:ir), kc + ( 1: ic) ] <- x [[i] ] 
kr <- kr + ir 
kc <- kc + ic 

} 

return (z) 

> 

repList <- function(x, n) { 
z <- list() 
for (i in l:n) 

z <- c(z, list(x)) 
return(z) 


shapeMe <- function (x) { 
m <- length (x) 

n <- (1 + sqrt (1 + 8 * m)) / 2 
d <- matrix (0, n, n) 
k <- 1 

for (i in 2:n) { 

for (j in 1: (i - 1)) { 

d [i, j] <- d[j, i] <- x [k] 
k <- k + 1 

> 

} 

return (d) 


symmetricFromTriangle <- function (x, lower = TRUE, diagonal = TRUE) { 
k <- length (x) 
if (diagonal) 

n <- (sqrt (1+8* k) -1) / 2 
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else 

n <- (sqrt (1+8* k) +1) / 2 
if (n != as.integer (n)) 
stop ("input error") 
nn <- 1: n 

if (diagonal && lower) 
m <- outer (nn, nn, ">=") 
if (diagonal && (! lower)) 
m <- outer (nn, nn, "<=") 
if ((Idiagonal) && lower) 
m <- outer (nn, nn, ">") 
if ((Idiagonal) && (! lower)) 
m <- outer (nn, nn, "<") 
b <- matrix (0, n, n) 
b[m] <- x 
b <- b + t(b) 
if (diagonal) 

diag (b) <- diag(b) / 2 
return (b) 


triangleFromSymmetric <- function (x, lower = TRUE, diagonal = TRUE) { 
n <- ncol (x) 
nn <- 1: n 

if (diagonal && lower) 
m <- outer (nn, nn, ">=") 
if (diagonal && (! lower)) 
m <- outer (nn, nn, "<=") 
if ((Idiagonal) && lower) 
m <- outer (nn, nn, ">") 
if ((Idiagonal) && (I lower)) 
m <- outer (nn, nn, "<") 
return (x [m]) 


5.3 mdsUtils.R 

library (mgcv) 

torgerson <- function (delta, p = 2) { 

z <- slanczos(-doubleCenter((delta 2) / 2) , p) 
w <- matrix (0, p, p) 
v <- pmax(z$values, 0) 
diag (w) <- sqrt (v) 
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return (z$vectors w) 


> 

basisPrep <- function (n, p, w) { 
m <- n * (n - 1) / 2 

v <- -symmetricFromTriangle (w, diagonal = FALSE) 

diag (v) <- -rowSums(v) 

ev <- eigen (v) 

eval <- ev$values[1:(n - 1)] 

evec <- ev$vectors[, 1: (n - 1)] 

z <- evec •/,*•/, diag (1 / sqrt (eval)) 

a <- array (0, c(n - 1, n - 1, m)) 

k <- 1 

for (j in 1: (n-1) ) { 
for (i in (j+1):n) { 
dif <- z [i,] - z[j ,] 
a [, , k] <- outer (dif, dif) 
k <- k + 1 

> 

} 

return (list (z = z, a = a)) 

} 

center <- function (x) { 

return (apply (x, 2, function (z) z - mean (z))) 

} 

doubleCenter <- function (x) { 
n <- nrow (x) 
j <- diag(n) - (1 / n) 
return (j */,**/, x j) 

} 

squareDist <- function (x) { 
d <- diag (x) 

return (outer (d, d, "+") - 2 * x) 


5.4 smacof.R 

smacof <- 

function (delta, 

P = 2, 

xini = torgerson (symmetricFromTriangle (delta, diagonal = FALSE), p), 
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w = rep (1, length (delta)), 
itmax = 1000, 
eps = le-10, 
verbose = FALSE) { 
m <- length (delta) 
n <- (1 + sqrt (1 + 8 * m)) / 2 
r <- p * (n - 1) 
h <- basisPrep (n, p, w) 
xold <- rep (0, r) 
for (s in l:p) { 

k <- (s - 1) * (n - 1) + 1: (n - 1) 
xold[k] <- 

crossprod (h$z, xini[, s]) / diag (crossprod (h$z)) 

> 

d <- rep (0, m) 
itel <- 1 
sold <- Inf 

ssq <- sum (w * delta ~ 2) 
repeat { 

bmat <- matrix (0, n-1, n-1) 
xx <- matrix (xold, n - 1, p) 
for (k in l:m) { 
ak <- h$a[, , k] 

d[k] <- sqrt (sum (xx * (ak xx))) 

bmat <- bmat + w[k] * (delta[k] / d[k]) * ak 

> 

xnew <- as.vector (bmat xx) 
snew <- sum (w * (delta - d) ~ 2) / ssq 
if (verbose) 
cat ( 

"Iteration: ", 

formate (itel, width = 3, format = "d"), 

"sold: ", 
formate ( 
sold, 

digits = 8, 
width = 12, 
format = "f" 

), 

"snew: ", 
formate ( 
snew, 

digits = 8, 
width = 12, 
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f ormat 


), 

"\n" 

) 

if ((itel == itmax) I I ((sold - snew) < eps)) 
break 

xold <- drop (xnew) 
sold <- snew 
itel <- itel + 1 

> 

xconf <- matrix (0, n, p) 
for (s in 1: p) { 

k <- (s - 1) * (n - 1) + 1: (n - 1) 
xconf [, s] <- h$z y,*y, xnew [k] 

> 

return (list ( 
delta = delta, 
dist = d, 
x = xconf, 
itel = itel, 
stress = snew 

)) 

> 


5.5 smacofBelow.R 

smacofBelow <- 
function (delta, 

P = 2, 

xini = torgerson (symmetricFromTriangle (delta, diagonal = FALSE), p), 
w = rep (1, length (delta)), 
bnd = delta, 
itmax = 1000, 
eps = le-10, 
verbose = FALSE) { 
m <- length (delta) 
n <- (1 + sqrt (1 + 8 * m)) / 2 
wnd <- which (bnd < Inf) 

1 <- length (wnd) 
r <- p * (n - 1) 
yold <- rep (0, 1) 
h <- basisPrep (n, p, w) 
xold <- rep (0, r) 
for (s in l:p) { 
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k <- (s - 1) * (n - 1) + 1: (n - 1) 
xold[k] <- 

crossprod (h$z, xini[, s]) / diag (crossprod (h$z)) 

> 

d <- rep (0, m) 
itel <- 1 
sold <- Inf 

ssq <- sum (w * delta ~ 2) 
repeat { 

bmat <- matrix (0, n-1, n-1) 
xx <- matrix (xold, n - 1, p) 
for (k in l:m) { 
ak <- h$a[, , k] 

d[k] <- sqrt (sum (xx * (ak °/ 0 *°/ 0 xx))) 
bmat <- bmat + w[k] * (delta[k] / d[k]) * ak 

> 

xmid <- as.vector (bmat °/ 0 *% xx) 

ccomp <- c(ssq,-bnd[wnd] 2) / 2 

bcomp <- matrix (0, r, 1 + 1) 

bcomp[, 1] <- -xmid 

acomp <- array (0, c(r, r, 1 + 1)) 

acomp[, , 1] <- diag (r) 

t <- 1 

for (k in wnd) { 

ak <- directSum (repList (h$a[, , k], p)) 
acomp[, , t + 1] <- ak 
t <- t + 1 

} 

v <- qpqc (yold, acomp, bcomp, ccomp, verbose = FALSE) 
snew <- 2 * v$f / ssq 
xnew <- v$xmin 
ynew <- v$multipliers 
cons <- v$constraints 
if (verbose) 
cat ( 

"Iteration: ", 

formate (itel, width = 3, format = "d"), 

"sold: ", 
formate ( 
sold, 

digits = 8, 
width = 12, 
format = "f" 

), 
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snew: 


formate ( 
snew, 

digits = 8, 
width = 12, 
format = "f" 

), 

"\n" 

) 

if ((itel == itmax) || ((sold - snew) < eps)) 
break 

xold <- xnew 
sold <- snew 
yold <- ynew 
itel <- itel + 1 

> 

xconf <- matrix (0, n, p) 
for (s in l:p) { 

k <- (s - 1) * (n - 1) + 1: (n - 1) 
xconf [, s] <- h$z 7.*7. xnew[k] 

> 

return ( 
list ( 

delta = delta, 
dist = d, 
x = xconf, 
multipliers = ynew, 
constraints = cons, 
itel = itel, 

stress = sum (w * (delta - d) ~ 2) / ssq 

) 

) 

} 
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